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Abstract 

The volume penalization method offers an efficient way to numerically sim- 
ulate flows around complex-shaped bodies which move and/or deform in 
general. In this method a penalization term which has permeability 77 and 
a mask function is added to a governing equation as a forcing term in or- 
der to impose different dynamics in solid and fluid regions. In this paper 
we investigate the accuracy of the volume penalization method in detail. 
We choose the one-dimensional Burgers' equation as a governing equation 
since it enables us extensive study and it has a nonlinear term similar to the 
Navier-Stokes equations. It is confirmed that the error which consists of the 
discretization/truncation error, the penalization error, the round-off error, 
and others has the same features as those in previous results when we use 
the standard definition of the mask function. As the number of grid points 
increases, the error converges to a non-zero constant which is equal to the 
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penalization error. We propose a new approach for reducing the penaliza- 
tion error, based on the analytical solution of the one-dimensional diffusion 
equation with a penalization term. We modify the mask function by shifting 
the boundary between the solid and fluid regions by ^JTn\ toward the fluid 
region. The numerical results for the one-dimensional Burgers' equation and 
the two-dimensional Navier-Stokes equations show that our new approach is 
effective in reducing the error; the error decreases to zero as the number of 
grid points increases. We also found that the error reduction is enhanced for 
large permeability. 

Keywords: volume penalization method, immersed boundary method, 
compact scheme, error reduction 



1. Introduction 

Flows around solid bodies have been studied in a wide variety of fields in 
science and engineering. In experiments it is difficult to visualize flow field 
and to obtain the field data in a three-dimensional space. In this regard com- 
putational fluid dynamics has advantages since the precise data are at hand. 
There are several kinds of space discretization methods, for example, the fi- 
nite difference method, the finite volume method, the spectral method, and 
so on. In the finite difference and the finite volume methods, a grid system 
is normally generated so that the boundaries between fluid and solid regions 
coincide with grid lines forming a body-fitted grid system because it is easy 
to impose boundary conditions. On the other hand, in the spectral method a 
set of orthogonal functions such as Fourier series and Chebyshev polynomials 
are used to expand the flow variables. It is used for flows in simple geome- 
tries like periodic boundary. However, if there exist complex-shaped solid 
bodies or bodies which move or deform in the flow, it is difficult to generate 
a body-fitted grid system of high quality and to simulate the flow using the 
spectral method since it is hopeless to find a set of orthogonal functions with 
computational efficiency The volume penalization (VP) method is one of 
the simplest methods that enable us to calculate flows which include solid 
bodies with mobility, deformability, and complex geometry. 

The VP method is one of the immersed boundary (IB) methods [ill ]. In 
the IB methods an external force term is added to the equations of motion 
to impose boundary conditions. IB methods are classified into two types: 
the continuous forcing approach in which an external force term is added 
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to a continuous equation and the discrete forcing approach in which it is 
added to a discretized equation. The VP method is the former type. We 
can use it with the Fourier pseudo-spectral method; many flows in which 
multiple solid bodies exist 0, 0, Qjl, 13 , thejlows inside rigid boundaries 



14|, |15|], and the flows around moving bodies [9| have been simulated by the 
VP method. Moreover, the VP method can be used with Chebyshev pseudo- 
spectral method, wavelet solvers, and other high-precision methods[6]. 

In the VP method, a solid body is regarded as porous medium. There 
are two types of penalization modeling. One is the L 2 penalization: the N-S 
equation is converted to the Darcy equation in the solid body; and the other 
is the H 2 penalization: the N-S equation is transformed to the Brinkman 
equation in the solid body [H, 0]. Most of the studies by VP methods have 
adopted the L 2 penalization. The added term, called a penalization term, 
acts as a damping force term which has a mask function \ and the perme- 
ability rj. Usually the step function, which is in the fluid region and 1 in 
the solid region, is chosen as \. The mask function activates the penalization 
term in the solid region, and the penalized N-S equation turns into the Darcy 
equation. 

One of the advantages of the VP method is that there are rigorous re- 
sults about convergence: the penalized solution converges to the solution of 
the original problem with no-slip boundary conditions as permeability tends 
to zero. Angot et a/.[l[ proved mathematically that the upper bound for 
the penalization error, which is the difference between the solutions of the 
original and penalized N-S equations, is 0(^ 1//4 ) in the fluid region. Carbou 
and Fabrie[i| showed that this upper bound is improved to 0(?7 1//2 ). More- 
over, Kevlahan and Ghidagliaj3] considered a Stokes flow over a flat plate 
whose dynamics is reduced to the one-dimensional diffusion equation and 
showed analytically that the error between the exact and penalized solutions 
is 0(r] 1 ^ 2 ) in the fluid region. These results imply that we can make the pe- 
nalization error smaller than the discretization error by choosing sufficiently 
small t]. 

However, we cannot choose too small value for rj when we use an explicit 
method for time development. The relation At < Crj should be fulfilled in 
order to ensure numerical stability, where At is the time step and C is a 
constant which depends on the method of time integration. When the total 
error, which is the difference between a penalized numerical solution and 
the exact solution to the original problem, is dominated by the penalization 
error, we cannot make the total error smaller than the penalization error in 
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the limit of large N which is the number of grid points. Keetels et a/.|6|] 
found that there is an optimal choice of N and rj giving the minimum total 
error when the penalization error approximately compensates the truncation 
error. They obtained the relation between the two parameters N and rj which 
give the optimal choice. If we know the two parameters in advance, we can 
make both the total error and the computational cost small by specifying 
the optimal choice of N and rj. However, their results show that the choice 
depends on the space discretization method. In addition, there is usually 
another requirement for the size of grid spacings that comes from the smallest 
scale of fluid motion. Therefore, other strategies are needed to reduce the 
total error when we cannot choose small value for rj. 

In this paper we investigate the accuracy of the volume penalization 
method. The characteristics of the error of the VP method are investigated 
in detail for the one-dimensional (ID) Burgers' equation using various space 
discretization schemes. Since exact solutions are known for the ID Burgers' 
equation in spite of its nonlinearity, we can evaluate the error between nu- 
merical and exact solutions explicitly. Next we propose a new approach for 
reducing the total error by modifying the mask function. Then we show that 
it is also effective for the two-dimensional (2D) N-S equations. 

This paper is organized as follows. First, in Section 2 we show the results 
of basic studies on the VP method and propose a new method for error 
reduction. Next, we apply it to a 2D problem in Section 3. Finally, we 
conclude in Section 4. 



2. Basic study on ID problems 

2.1. Statement of the problem 

In this section, we are mostly concerned with the ID Burgers' equation 

du du d 2 u 
dt dx dx 2 ' 

with the initial and boundary conditions 

u(0, x) = u (x) = - sin (J^-J , (2) 

u(t,±L) = 0, (3) 

where v denotes the diffusion coefficient. The function u(t, x) may be re- 
garded as flow velocity as Eq. ([I]) is the ID compressible N-S equation without 
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pressure. Kevlahan and Ghidaglia[3] considered the ID diffusion equation in 
order to investigate basic properties of the VP method. The ID diffusion 
equation is desirable for validating numerical methods since it is linear so 
that both the original equation and the penalized equation have analytical 
solutions. However, the effects of nonlinearity cannot be tested. Thus we 
choose the simplest nonlinear equation which is related with the N-S equa- 
tions and has analytical solutions for the original equation. 

Exact solutions of Eqs. ([I])-© are obtained using the formula by Cole 

u t (t x) = 2i/ Sijkl!^Z^ exp[-un 2 n 2 t/L 2 }A n sm(nnx/L) 
A + Y^Li ex P[ — vn 2 n 2 t/L 2 ]A n cos(nnx/L) 

A) = 777 / 9o(x)dx, (5) 



2L 



-L 



1 f L „ . , n-nx 



A n = — J 6 (x) cos— j—dx, (6) 

9 (x) = Cexp (-l.J\(Z)dtJ . (7) 

Using the VP method, we try to obtain a numerical solution as close as 
possible to the analytical solution (j3J). We set L = tt in the rest of the paper. 

2.2. Numerical setups 

In the VP method, the ID Burgers' equation in Eq. (pQ) with the boundary 
conditions (E} is approximated by the equation with a penalization term 

du du d 2 u X ( \ ( o\ 
dt dx dx 2 7] s ' 

where u s is prescribed velocity of a solid region and 77 is permeability. In this 
paper, u s is basically zero, v is fixed to be 10 _1 , and rj is varied from 10~ 2 to 
10 -5 . The mask function \ is se f f° a sf e P function, 

^ { : : ;:: ::: ^ 

The whole computational domain Q = Qf + Q s covers the range |x| < Lf,. 
The fluid region flf is defined as \x\ < L, and fl s corresponds to the solid 
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regions L < \x\ < L^, where Lf, = 2tt. Thus the boundaries between Qf and 
fl s are located at x = ±L = ±7r. The initial condition is 

( \ / -sin (a:) in Q f 
«(0,z) = | in Q [ ■ (10) 

The terms containing spatial derivatives in Eq. (jSJ) are discretized with 
one of the following methods: the second-order central difference method, the 
fourth-order central difference method, the fourth-order Pade type compact 



finite difference scheme [10J, and the Fourier spectral method. We use a 
uniform mesh. The fourth-order Runge-Kutta method is used for the time 
marching. The time step is fixed to be At = 10~ 5 . The grid/mode number of 
N = 250 — 2750 and the permeability of rj = 10 -5 — 10~ 2 are chosen to meet 
the conditions due to the nonlinear term u maK At/Ax < C c , the diffusion 
term uAt/(Ax) 2 < C4, and the penalization term At < Crj, where M max is 
the maximum velocity and C c , Cd, and C are constants. 

Here we consider the definition of the step function on discrete grid points. 
We cannot express the step function rigorously on discrete grid points since 
there is a non-zero gap between the grid points at which the value jumps from 
to 1. Thus we should choose a mask function which gives the penalized 
numerical solution correctly converging to the penalized exact solution as 
iV tends to 00. Fig. [1] depicts three candidates for mask functions near the 
boundary between Q s and Qf\ the boundary coincides with the grid point 
where x is 1 011 Type A (Fig. [T^i); the boundary is located at the midpoint 
between the two grid points at which the value of the mask function jumps 
from to 1 on Type B (Fig. [T] b); the boundary coincides with the grid 
point where x is on Type C (Fig. [lb). In Section 1231 we first compare the 
results on Type A with previous results since it is a natural choice. Then, we 
determine which type of the definition of the mask functions is most suitable, 
and show how the errors are reduced when we modify the mask functions 
defined by the most suitable definition in Section 12.41 

2.3. Error analysis 

In this section, we investigate the characteristics of the total error for 
Type A mask function. After showing the general features of the total error 
of the penalization method, we compare them with previous results. 

2.3.1. General features 

Fig. [2] compares the numerical and exact solutions of the ID Burgers' 
equation in the whole computational region (Fig. |2h) and near the boundary 
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between Q s and Qf (Fig. [2b). The numerical solutions are computed with 
the highest resolution until t = 100. The solutions are close to a sine wave in 
the fluid region and is approximately zero in the solid regions (Fig. |2ji). For 
large values of rj, however, the numerical solutions deviate from the exact 
solution (Fig. [2^,) and there is leak from solid to fluid regions if we regard u 
as flow velocity (Fig. [2b). This deviation from the exact solution is due to 
the penalization term. They decrease as rj tends to zero. 

Next, we discuss the grid/mode number dependence of total errors using 
(a) the second-order central difference method, (b) the fourth-order central 
difference method, (c) the fourth-order compact difference method, and (d) 
the spectral method. In this study, total error 5 tot is defined as the root mean 
square of the difference between a numerical solution and the exact solution 



J Qf \u(t,x) - u cxact {t,x)\ 2 dx 
In f ^ x 

where w ex act is an exact solution of the ID Burgers' equation without the 
penalization term, obtained from Eqs.(j4]) - (|7|). The total error <5 tot may be 
expressed as 

5tot = 5 V + 5 N + 5 etc , (12) 

where 5 V , 5^, and 5 etc are the penalization error, the space discretization/truncation 
error, and the sum of the other errors. Note that 5 V and 5^ are independent 
of iV and 77, respectively, and 5 etc is normally negligible. The data plotted 
by asterisks in Fig. [3] are obtained by solving the non-penalized ID Burgers' 
equation with the no-slip boundary conditions. As shown in Fig. [3] (d) for the 
spectral method 6 to t takes a constant value of O(10~ 14 ), which is regarded as 
the value of 5 etc - 

In Figs. [3J (b), (c), and (d), the optimum grid/mode number N opt at 
which 5 tot takes the minimum value is clearly observed for each value of 77. 
At N = N opt , <5jv and 5 V would have the same magnitudes and opposite signs 
to be canceled out. For N > N opt the total error S tot tends to a constant with 
increasing N. In other words <5 tot does not depend on N when iV 3> iV opt , 
implying that 5 tot is dominated by 5 V . For iV A^ opt the total error for all 
values of 77 converges to one line, implying that <5 tot does not depend on 77 for 
iV iV opt , as <5 to t is dominated by 5^. The dependence of S to t on N shown 
above agrees with the previous results by Keetels et al. 0]. In the following 
three sections, we investigate the features of these errors in more detail. 
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2.3.2. Optimum grid/mode number at N = N opt 

Keetels et al. jsf performed an error analysis of the VP method for the 
2D N-S equations by numerical simulation of the dipole-wall collision using 
the pseudo-spectral method. They found relations 5n = o>N a and 5 V = (3r] b , 
where a, (3, a(< 0), and 6(> 0) are constants which depend on numerical 
schemes. Previous works show b = 0.50, @]. In addition, the following 
relation is obtained 



where 7 and c are also parameters, by setting 5^ = — 5 V which makes 8 tot 
small. For a given value of 77, the value of N determined by Eq. f|T3|) corre- 
sponds to iV p t . In Figs. [3] (b), (c), and (d), iV p t increases with decreasing 
T], which is consistent with Eq. ( {TBI since c = a/b < 0. Keetels et al. |6[ 
suggested that we can perform the most accurate and efficient simulation if 
we can determine the combination of N and 77 by Eq. ({TBI . To do this, we 
must specify the unknown constants 7 and c. 

Here, we obtain 7 and c for each discretization method by substituting 
iV pt into N in Eq. ({TBI) . The relation between 77 and iV pt is shown in Fig. HI 
The higher-order schemes have the same power c = —2.0, but they have 
small differences in 7. As seen in Fig. [3], the ranges of N giving small error 
are too narrow to allow the small differences between the values of 7 for 
high-order schemes. For instance, if we obtain iV p t from the results of the 
spectral method and use it for simulations with a different discretization 
scheme, the resulting error 5 to t is not small. Moreover, 7 and 77 for the 
second-order central scheme are rather different from those for the higher- 
order schemes. We infer from Fig. [3] (a) that the second-order central scheme 
cannot accomplish Sn = — 8 V . 

Next, we investigate the dependence of 7 and c on the problem. The 
values of 7 and c in the present and previous studies are summarized in the 
Table [TJ The left column shows the present study. The middle column shows 
the results of the parallel plate channel flow governed by the ID diffusion 
equation^]. The right column shows the results of 2D simulations of dipole- 
wall collision by solving the N-S equations by the spectral method or the 
coherent vortex simulation (CVS) [6]. Comparing the values for the second 
central scheme or the spectral method, it is clear that 7 and c depend not 
only on the spatial discretization but also on the problems. 

The above results show that it is difficult to decrease 5 to t using Eq. (113"]) 




(13) 
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which includes unknown parameters 7 and c, because they are so sensitive to 
both the space discretization scheme and the problem that we cannot specify 
them to get a benefit of the error reduction. Thus we need another strategy 
for decreasing 5 to t, which will be proposed in Section 12^1 

2. 3. 3. Penalization error for N opt <C N 

The penalization error dominates the total error for N 3> N opt . Fig. [5] 
(a) shows the total error 5 tot in the fluid region Qf as a function of r] for 
the fourth-order compact scheme. The total errors are the same ones for 
N m3iX = 2500 shown in Fig. [3] (c), which are chosen as a typical example. 
The total error 5 to t tends to be proportional to 77 ' 5 as rj increases, where N 
is sufficiently larger than iV p t . The characteristics agrees with the analytical 
results for the dependence of the penalization error on r\ for the N-S equation 
[i| and the diffusion equation Q. 

We also show <5 tot as a function of 77 in the solid region fl s (Fig. [5b). It 
is calculated by Eq. ffTTl) replacing Qf by Q s . In the solid region, 5 to t does 
not depend on iV very much, although there are some differences between 
the results for the minimum grid number iV min = 300 and for the maximum 
iV max = 2500. The dependence of 5 tot on r\ varies from 0(7? a82 ) to O(?7 ' 98 ). 
The same results are obtained numerically for a flow past a square cylinder 

As for the other space discretization methods, the results similar to above 
are obtained for the penalization error. To summarize the present results 
obtained for a nonlinear equation show the same features as those obtained by 
the theoretical studies [i| and the numerical studies for both the ID diffusion 
equation 0] and the multi-dimensional N-S equations However, as far as 
we use an explicit method for time evolution with a time step which is not 
sufficiently small, the penalization error is much larger than the error of the 
conventional method of imposing boundary conditions on boundary-fitted 
grids (Fig. [3]). This is another reason why we seek a strategy for reducing 
the penalization error. 

2.3.4- Discretization error for N <C iV pt 

As can be seen in Figs. [3] (b), (c), and (d), the total error converges to a 
single asymptotic line which is independent of 77 as N decreases. This implies 
that the discretization error is the dominant component of the total error for 
iV <C iV p t . The asymptotic line is proportional to iV -1 showing that the 
spatial accuracy of the simulation is first order. The convergence of the total 
error for the ID channel flow by Kevlahan and Ghidaglia [7| also indicates 
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the first-order accuracy unless N rs N opt . According to them this low-order 
accuracy is due to the C 2 discontinuity of the penalized exact solution at the 
boundaries. 

In order to check it, we consider a steady solution of the ID Burgers' 
equation 

u(oo,x) = — Mbtanh [t^x^i (14) 

which satisfies u(t, x) — > T^b, du(t,x)/dx —> 0, and d 2 u(t,x)/dx 2 — > for 
1 <C Ub\x\/2u. For this solution the effect of penalization can be minimized 
if the boundaries between fluid and solid regions are located sufficiently far 
from x = 0. We solve Eq. (JSJ), where we set 

{Ub if x < 
if x = . (15) 
— Ub if x > 

The mask function (Q is unchanged. The initial condition is given by the 
steady solution f lT4|) . 

/rt >. f — it& tanh ( — x ) in Qt 
m(0,x) = < V2z/ / f , (16) 

[ — sgn(x) in fi s 

where = 1. We use the numerical solutions which have converged with 
sufficient accuracy for the error analysis below. 

Fig. [6] shows the steady solutions of the ID Burgers' equation numeri- 
cally obtained using the fourth-order compact scheme with iV = iV max . The 
penalized numerical solutions for rj = 10~ 2 and 10~ 3 agree well with both 
the non-penalized numerical solutions and the exact solutions (Fig. E^). The 
gradients of the exact solution at the boundaries are almost zero (Fig. |Bb). 
That is, the error due to the discontinuity at the boundaries is within 5 etc 
in Eq. ([P2]) . Thus the numerical solutions can be regarded as sufficiently 
smooth in this problem. 

Fig. [7] shows the total error as a function of iV for the steady problem. 
Both the total error 5 to t for the penalized numerical solution and that for the 
non-penalized numerical solution with boundary conditions decrease as iV -4 . 
The accuracy of the spatial discretization scheme used in this simulation is 
achieved when the solution is smooth at the boundaries. Thus we conclude 
that the discontinuity of the exact penalized solution is responsible for the 
low-order accuracy of the discretization error of the penalization method. 
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The accuracy of the VP method may be improved by using a continuous 
mask function instead of the step function. However, it is not clear whether 
the solutions of the penalized equation with a continuous mask function 
converge to that of the non-penalized equation as rj tends to zero. In this 
paper, we focus on the reduction of the penalization error while the reduction 
of the discretization error is not pursued. 

2-4- Error reduction method 

In this section, we propose a new method for reducing the penalization 
error by changing the mask function. First, we determine the appropriate 
definition of the mask function for good convergence. Next, we modify the 
mask function so that the exact solutions of the penalized equation almost 
agrees with that of the non-penalized equations in the fluid region. 

2.4-1- Most suitable definition for mask functions 

We have introduced three types of mask functions: the values at the 
boundaries are x(±L) = 1.0 for Type A, 0.5 for Type B, and 0.0 for Type 
C (Fig. [T]). The total error as a function of grid/mode number for the three 
types is shown in Fig. [HJ We found that the values of 5 to t are different among 
the three types. For Types B and C we note that there is no iV pt at which 
5 to t takes local minimum value in the region 300 < N < 2750. That is, the 
numerical solutions for Type B and C cannot satisfy any conditions to meet 
8n — — <V The above observation suggests that Keetels et al. [6| used the 
mask functions of Type A since iV p t existed in their results. However, there 
are also common features among the three types. For small N, the total error 
<5 tot increases and collapses to a single line with decreasing N regardless of 
rj. For large N, on the other hand, 5 tot converges to a value which depends 
on rj as iV increases. The latter fact implies that the penalized numerical 
solution does not converge to the non-penalized exact solution w exact . Instead 
it converges to the penalized exact solution u v ex act- Therefore, we should 
consider the difference between the numerical solution and u v cxa ct in order 
to compare the convergence properties among the mask functions for the 
three types. 

Since it is difficult to obtain u v ex act, we introduce 




L \u(t,x) ~u Nrn ^(t,x)\ 2 dx 



(17) 
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where u^max is the numerical solution for N = iV max , to examine the conver- 
gence property of numerical solutions. If N max 3> N, «Ar max can be regarded 
to be close to u v ex act- Fig. [9] shows the error (5jv max defined above as a function 
of N. For all three types of mask functions Spf max decreases with increasing N. 
Therefore, regardless of the type of a mask function, the penalized numeri- 
cal solution monotonically converges to some solution. For iV <C iV max the 
error 5jv max is proportional to iV~ L0 for Type A, A^~ 2,0 ~ A^~ L0 for Type B, 
and A^ -1 ' 3 for Type C. For Type B we have second-order accuracy for large 
T], which is the highest accuracy among the cases considered. Furthermore, 
we notice that as 77 increases decreases for Types B and C, while it 

increases for Type A. 

To clarify the above results, we calculate the difference between the nu- 
merical solution and the penalized exact solution u v ex act for a simpler prob- 
lem. Let us consider the penalized ID diffusion equation 

du = ^(Pu _ x u , lg v 
dt dx 2 rj ' 

instead of the penalized ID Burgers ' equation, whose exact solution is diffi- 
cult to find analytically. The exact solution of Eq. f ll8p is 

C+e -uk' n h e -ax L < X 

u v exact (t,x) = <j C e~ uk '* 1 sin (k' n x) \x\ < L , (19) 

where k' n is the wavenumber of a solution for the penalized ID diffusion 
equation and C± and a are determined by the condition of C l continuity 



v exact (t, ±L ~ 0) = U n exact ±L + 0), (20) 



U 

^y xact (t, ±L - 0) = ^ cxact (t, ± L + y (2i) 

By setting Co = — 1, the constants are found to be 

C± = Te aL sm{k' n L), (22) 

a_ = cos (k' n L) 
K sin (k' n L) ' 
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In addition, on substituting Eq. ( TIT)]) into Eq. (ITS]) , we obtain 



(24) 

We numerically solve the set of equations (123]) and (|2"%]I to obtain a and 
fc^j. We assume k' n be 0(1) while L is set to it. Substituting a and k' n into 
Eq. (Tl9]) . the solution exact in Eq. ( [18]) is determined. 

Next we define 5 V ex act as the difference between the numerical solution 
and the penalized exact solution 



"fj exact 



• r ) — M ^ exact(^) x)\ 2 d,X 

\1 ]^T X ' (25) 



The profiles of 5 V cxact as a function of AT are depicted in Fig. [TO] When 
iV is small, 5 n ex act for the penalized ID diffusion equation has the same 
convergence properties as for the penalized ID Burgers' equation as 

shown in Fig. [9] Thus, we can estimate how 5 V exa ct for the penalized ID 
Burgers' equation converges in the large limit of N by observing 5 n cxact 
for the penalized ID diffusion equation. For Type B 6 V exact tends to be 
proportional to iV -2 for large N, although it is proportional to N^ 1 for 
small N. On the other hand, for Types A and C, the decay rate of 5 V exact 
decreases for large N. The scaling of 5 rj cxa ,ct for Types A and C is flatter 
than 0(A^~ 1 ) and 0(A^~ 12 ), respectively. Therefore, the best definition of 
the mask function is Type B, which makes the numerical solutions converge 
to u v exact with second-order accuracy. In the following, we discuss in detail 
the accuracy and the dependence of 5 V exact on 77 using the mask functions 
Type B. 

The exact solution of the penalized ID diffusion equation has the C l 
continuity because of Eqs. ( 120]) and f ]2T]) . However, the second derivative at 
the boundaries between fi/ and fl s is discontinuous 

d 2 u , , , , d 2 u 



()l J^ ±L -°)^^ ±L + °)- ( 26 ) 

The ID diffusion equation (ITS]) has a second derivative term. This is why 
the accuracy is second-order although fourth-order schemes are used for the 
spatial discretization. The solutions of the penalized ID Burgers' equation 
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and the N-S equations are also inferred to be same as those of the penalized 
ID diffusion equation. 

Fig. [TT1 shows the solutions near the boundary for Type B. The solid line 
with plus symbols is the non-penalized exact solution u exa ct- The solid lines 
with other symbols denote the penalized numerical solutions u and the other 
lines are the penalized exact solutions u v exac t- Both u and u n exact approach 
""exact as r] decreases. However, u deviates from u v exact as r\ decreases. The 
penalized solution has a leaching area in a solid region. The width of the 
leaching area is 

l/a « v^tj, (27) 

which is estimated by substituting Eq. ( |24l) into Eq. ( Tl9l under the condition 
^Jvfj <C 1. Thus the leaching area becomes smaller with decreasing 77 as can 
be seen in Fig. [TTJ The leaching area for small 77 is too small to be resolved 
with N = iV max = 2750 grid points. It is the reason why 5 V exac t is large for 
small 77 (Fig. ITOb). 

2.4-2. Modified mask functions 

Substituting a and k' n obtained above into Eq. (fl9j) . the exact solution of 
the penalized diffusion equation is expressed as 



e -uk n t gin Q^ L j e -a(x+L) L < 



X 



w„exact(*,a;) = { -e- vk ^sm{k' n x) \x\ < L . (28) 

+e -K 2 t S i n [k' n L) e +a( - x+L ^ x < -L 

The exact solution of the non-penalized diffusion equation is 

f L < x 

u exac t(t,x) = < -e - "^ sin (k n x) \x\ < L , (29) 
[0 x < —L 

where k n — ^ (n = 1, 2, • • ■). The wavenumber fc^ in Eq. (l2~8~i) corresponds to 
the closest k n . We set L = n. Comparing Eqs. (1281) and (1291) . it is legitimate 
to shift the mask function replacing L by Lyp so that k' n = k n . Then the 
relation (123]) turns out to be 

sin k n Ly P k n 
cos k n Lyp a 

Substituting 



L + e, (31) 
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into Eq. (I50|) we obtain 

sin (k n L) cos (k n e) + cos (k n L) sin (k n e) k n 
cos (k n L) cos (fc n e) — sin (k n L) sin (/c n e) a 

Assuming \k n e\ <C 1 Eqs. (I2TJ) and ( 1321) give 



(32) 



e « -y^- (33) 

Therefore, by setting the boundaries at x = ±Lyp = ±(L — we obtain 
the penalized numerical solutions of which wavenumber k' n almost coincides 
with k n . 

The original and shifted mask functions for Type B near the solid bound- 
ary are shown in Fig. [12j The original boundaries are located at x — ±L = 
±7r, shown by the thin broken lines. The shifted function is 

x<») = { ; % . (34) 

where Q'f is \x\ < Lyp and Q' s is Lyp < \x\ < Lb. The shifted boundary 
x = —Lyp is shown by the thin solid line in Fig. [12] (b). Since the distance 
between the grid points is adjusted to shifted the mask function, the overall 
numerical domain Q', which is |x| < L&, depends on N although is close 
to 2tt. 

The total error S' tot for the shifted mask function is defined by Eq. (lllj) . 
replacing Qf by Q'f. The total error S' tot is shown as a function of N' for 
the penalized ID diffusion equation in Fig. [131 where N' = L/^L^/N) is the 
grid/mode number. The error 5' tot is almost same with 5 V exa ct- By modi- 
fying the mask function the total error decreases as iV increases instead of 
converging to a non-zero constant; that is, the penalized numerical solutions 
converge to the non-penalized exact solutions. Since the ID diffusion equa- 
tion is linear the Fourier modes of the solutions evolve independently of each 
other. Thus we can simply regard that each mode whose wavenumber is k n 
is converted to the corresponding mode whose wavenumber is k' n by shifting 
the boundaries between fluid and solid regions by y/vrj. 

We verify whether the above shifted mask function is also effective for 
the nonlinear ID Burgers' equation. Fig. [TH shows the total error S' tot as 
a function of N'. The convergence of 5' tQt for the ID Burgers' equation 
has the same characteristics with that for the diffusion equation shown in 
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Fig. [121 The penalized numerical solutions obtained by using the shifted 
mask functions converge to the exact solutions of the ID Burgers' equation 
for — Lyp < x < Lyp. The nonlinear term in the ID Burgers' equation makes 
the Fourier modes interact with each other producing the high wavenumber 
modes from the low wavenumber modes. Nevertheless, the new technique 
of modifying the mask function gives the same effect on the numerical solu- 
tions of the nonlinear Burgers' equation. Compared with the case of normal 
mask functions in Fig. [8] (b), the total error for the shifted mask functions 
is reduced for large N, especially for large r\. A possible reason for the r\ 
dependence is discussed in the last paragraph of the section 12.4.11 

Finally, we discuss limitation of our approach. The solutions of the ID 
Burgers' equation are shown for various methods in Fig. [T5J The dashed lines 
are the original boundaries ±L between Q s and Qf, and the solid line is one of 
the shifted boundaries ±Lyp between Q r s and Q'*. In Fig. IToT a). the numerical 
solution for the shifted mask function agrees with the exact solution of the 
non-penalized ID Burgers' equation, while the numerical solution for the 
normal mask function deviates from the exact solution. Close look at the 
boundary shown in Fig. [TS] (b) reveals that the numerical solution for the 
shifted mask function agrees well with the exact solution for — Lyp < x < 
Lyp. However, there is a small difference between them in the vicinity of the 
boundaries for Lyp < \x\ < L. This difference cannot be eliminated by our 
new approach. 

3. Application to 2D problem 

In the previous section we have proposed a new method to reduce the pe- 
nalization error. In this section we apply it to the penalized 2D N-S equations 
and see whether it is effective. 

3.1. Numerical setups 

We perform numerical simulation by solving 2D N-S equations with a 
penalization term 

+ Vu = -— Vp + z/V 2 u- -(u-u s ), (35) 
ot po V 

V ■ u = 0, (36) 

where u = (u,v), and u s = (u s ,v s ), are the flow velocity and velocity in 
the solid region, respectively. In Eq. ( )35|) . p is pressure, po is density which 
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is constant, v is fixed to 10 2 , and r\ is set to 10 2 , 10 3 , and 10 4 . The 
normal-type mask function is 

*M = { \ Zni ■ < 37 > 

and the shifted-type mask function is described as 

*w (38) 

The polar coordinate system (r, 0) is also used in the following. We 
consider the Taylor-Couette flow between two co-axial cylinders. The inner 
and outer cylinders have the radii R\ and R2, and they rotate with the 
angular velocities u\ and U2, respectively. The overall computational domain 
Q is —71 < x,y < 7r. The fluid region Qf is defined as R\ < r < R2. The 
solid region Q s consists of the inner solid region fi sl = {(x,y) | < r < Ri} 
and the outer solid region Q s2 = {(x, y) \ R 2 < r E ft}. For the shifted- 
type mask function, we define Q'j = {(x,y) \ R\ + y/vrj < r < R 2 — y/vrj}, 
M'si = {(%,y) \ < r < Ri + y/vrj}, and fi' s2 = {(x, y) \ R 2 — ^/urj <r eCl}. 

The azimuthal component of velocity ue s is set to 

in n f 

u es {x) = { rui in Q sl , (39) 
rco 2 in fi s2 

and the radial component u rs , is zero all over the computational domain. 
The Taylor-Couette flow is given by 

U2R\-u x R\ {u x -U2)R\R\\ . 

2 2 — r ^ 2 2 111 ^/ 

u e (0,x)={ R2-R1 r (4Q) 

' ruji m S! sl 



rw 2 in Vt 



«2 



which is used as the initial condition, while the radial component u r (0,x) is 
zero in the whole computational domain. [16|. In this study, we set Ui = 1, 
U2 = 0, Ri = 0.4vr, and R 2 = 0.8vr. 

The Fourier pseudo-spectral method and the fourth-order Runge-Kutta 
method are used for spatial and time discretization, respectively. The mode 
number iV is varied from 256 to 4096. The time step At is 10~ 4 or 10 -5 
depending on N to meet the numerical stability condition. 
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Fig. [TBI shows the mask functions near the inner cylinder at various angles. 
Because the boundaries between fl s and flf (ft' s and Q'*) are circular, the 
radial position of the actual boundary depends on the angle. 

3.2. Error analysis and validation 

In this section, we examine the characteristics of the total error for the 
2D problem, and demonstrate the general applicability of our new approach. 

The total error <5 to t as a function of N for the normal and shifted mask 
functions are shown in Fig. [T71 The total error 5 to t calculated by Eq. ( TTTj) 
decreases as N increases showing second-order accuracy for the shifted mask 
functions, while it converges to a constant value for large N for the normal 
type mask functions. Moreover, the effect of error reduction increases with 
i]. It is the same feature as that for the ID diffusion and Burgers' equations 
(Figs. H~3l and [T4"j) . This result shows that our new method for error reduction 
is also effective for the 2D N-S problem. 

The solutions in the whole region, near the inner boundary, and near 
the outer boundary are shown in Fig. [TSJ The Taylor-Couette flow is well 
resolved by the pseudo-spectral method in all cases (Fig. IT8r). Looking 
at the vicinity of the inner and outer cylinders, however, we observe some 
differences between the numerical and the exact solutions (Figs. [THb and c). 
The results obtained by the normal type mask function are smaller /larger 
near the inner/outer boundary than the exact solution. By shifting the mask 
function toward the fluid region by y/vfj, the numerical solution agrees well 
with the exact solution in the fluid region except for the immediate vicinity of 
the cylinders. These features are exactly the same as those of the ID Burgers' 
equation. In addition, the radial distributions for 6 = 0° and 9 = 45° are 
in good agreement with each other for iV > 1000. Thus, the effect of error 
reduction is independent of the azimuthal angle for large iV although the grid 
system is not axisymmetric. 

4. Summary 

We have investigated the error of the volume penalization method, and 
have proposed a new method for reducing the error due to the penalization 
term. Our findings are summarized below. 

First, we analyzed the error of numerical solutions of the ID Burgers' 
equation with a penalization term. We used the mask function Type A, whose 
value at the boundaries is 1.0. The features of the total error as a function of 
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grid/mode number N are almost the same as those in previous works. There 
is an optimum grid/mode number A^ opt depending on the permeability r\ at 
which the total error takes a minimum value. The total error decreases as N 
increases for N <C iVopt, in which the discretization error is dominant. On the 
other hand, the total error becomes constant as N increases for N 3> N opt 
in which the penalization error is dominant. We tried to find how to specify 
iV p t to reduce the error. However, it turned out to be difficult to specify N opt 
in advance, because N opt depends on both the space discretization method 
and the problem under consideration. 

Next, we modified the mask function in order to reduce the total error. 
We found that the mask function Type B, for which the boundary is located 
at the midpoint of the grid points where the mask function jumps from to 
1, makes the numerical penalized solutions converge to the exact penalized 
solutions. Although Keetels et aZ.jfj] seem to have used the mask function 
Type A, it is important to use the mask function Type B to reduce the 
penalization error. Our new idea is to shift the boundary between fluid and 
solid regions of the mask function toward the fluid region by ^Jurj. The 
modified mask function makes the total error decrease as N increases so that 
the penalization error is significantly reduced. Since the leaking area in the 
solid region should be adequately resolved for given N to obtain the effect 
of error reduction, this technique is effective for comparatively large 77. Thus 
it is useful when the explicit method is used for time integration because of 
the condition for numerical stability At < Crj. 

Finally, we performed two-dimensional numerical simulations of Taylor- 
Couette flow between two co-axial cylinders to confirm the applicability of 
the present method. The results showed that modifying the mask function is 
also effective for two-dimensional incompressible Navier-Stokes equation even 
if circular solid boundaries are immersed in the Cartesian grids. Therefore, it 
would be valid for various governing equations, space discretization methods, 
and multi-dimensional problems. 

There are several conditions for permeability 77, some of which are men- 
tioned above: (i) time resolution: response time 77 should be smaller than the 
smallest time scale which should be resolved; (ii) spatial resolution: the "sur- 
face thickness" y/vrj should be smaller than the smallest length scale which 
should be resolved; (iii) resolution at the boundaries: the grid spacing should 
be smaller than s/vfj to resolve the surface "layer" at the boundaries; (iv) 
numerical stability for explicit time integration: At < Crj. The conditions 
(i) and (ii) are necessary, while (iii) is optional and (iv) is irrelevant when an 
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implicit method is used for time integration. 

In this study, we did not discuss the applicability of our error reduction 
method for the non-uniform grid, moving or deforming solid boundaries, vol- 
ume penalization method for compressible flow, and so on. Furthermore, we 
do not deal with continuous mask functions [9j]. These problems are important 
for development of the volume penalization method and will be investigated 
as future works. 
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Tabic 1: Summary of the values of 7 and c for the relational expression of 77 = jN£ pt in 
the present and previous studies. 





ID Burgers 
sine wave 


ID diffusion 
parallel plate channel 
by Kevlahan (2001) 


2D N-S 

dipole-wall collision 
by Keetels (2007) 


2nd central 


7 ~ 13000 
c = -3.1 


c = -3.0 




4th central 


7 ~ 110 
c = -2.0 






4th compact 


7 ~ 130 
c = -2.0 






Spectral 
Fourier 


7 ~ 155 
c = -2.0 


c = -2.3 


7 ~ O(10 5 ) 
c = -3.6 


CVS 
Wavelet 






7 ~ O(10 2 ) 
c = -1.5 
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Figure 1: The mask function for (a) Type A, (b) Type B, and (c) Type C. 
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Figure 2: The distribution of solutions for the ID Burgers' equation (a) in the whole 
region and (b) near the solid boundary. 
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Figure 3: The total error <5 t ot as a function of N for the various discretization methods: (a) 
the second-order central difference scheme; (b) the fourth-order central difference scheme; 
(c) the fourth-order compact difference scheme; (d) the spectral method. 
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Figure 4: The relation between the permeability i] and the optimum grid/mode number 
iV opt . 
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Figure 5: The total error 5tot as a function of 77 in (a) the fluid region J7/ and (b) the 
solid region fl s . 
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Figure 6: The steady solutions of the ID Burgers' equation (a) in the whole region and 
(b) near the solid boundary. The numerical solutions are obtained using the fourth-order 
compact scheme. 
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ID steady Burgers, 4th compact, TYPE A 
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Figure 7: The total error <5 to t as a function of N for the steady numerical solutions of ID 
Burgers' equation solved with the fourth-order compact scheme. 
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Figure 8: The total error #tot as a function of N for the numerical solutions of the ID 
Burgers' equation using the mask function for (a) Type A, (b) Type B, and (c) Type C. 
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Figure 9: The error <5/v max between the numerical solutions at N and iV max for the ID 
Burgers' equation using the mask functions for (a) Type A, (b) Type B, and (c) Type C. 
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Figure 10: The error 8^ exact between the numerical solutions u and the exact penalized 
solutions u v cxac t for the ID diffusion equation using the mask functions for (a) Type A, 
(b) Type B, and (c) Type C. 
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Figure 11: The numerical and exact solutions for the ID diffusion equation using 
Type B mask function. 
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Figure 12: The normal and shifted mask functions for Type B near the solid boundary. 
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Figure 13: The total error <^ ot as a function of N' for the ID penalized diffusion equation 
with shifted mask functions for Type B. 
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Figure 14: The total error 6' tot as a function of N' for the ID penalized Burgers' equation 
with the shifted mask functions for Type B. 
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Figure 15: The comparison of the solutions for the ID Burgers' equation (a) in the whole 
region and (b) near the solid boundary. 
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Figure 16: The distributions of x( x ) near the inner cylinder for the 2D problem: (a) the 
normal mask functions, (b) the shifted mask functions. 
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Figure 17: The characteristics of <5tot as a function of N for the 2D N-S equations: 
the normal mask functions, (b) the shifted mask functions. 
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Figure 18: The comparison of solutions for 2D N-S equations. The radial distributions 
of azimuthal velocity are shown (a) in the whole region, (b) near the inner cylinder, and 
(c) near the outer cylinder. 
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